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ABSTRACT 

We construct self-similar, axisymmetric, time-independent solutions to Ein- 
stein's field equations for an isothermal gas with a flat rotation curve in the 
equatorial plane. The metric scales as ds'^ a^ds"^ under the transformation 
r — s> ar and t a^~'^t, where n is a dimensionless measure of the strength of the 
gravitational field. The solution space forms a two-parameter family character- 
ized by the ratios of the isothermal sound speed and the equatorial rotation speed 
to the speed of light. The isodensity surfaces are toroids, empty of matter along 
the rotation axis. Unlike the Newtonian case, the velocity field is not constant 
on a cylindrical radius because of frame dragging. As the configuration rotates 
faster, an ergoregion develops in the form of the exterior of a cone centered about 
the rotation axis. The sequence of solutions terminates when frame dragging be- 
comes infinite and the ergocone closes onto the axis. The fluid velocity of the 
last solution has a modest value in the midplane but reaches the speed of light 
on the axis. 



1. Introduction 



In remarkable treatments of the axisymmetric equilibria of self-gravitating, isothermal, 
unbounded, stellar and gaseous systems with flat rotation curves, Toomre (1982) and Hayashi 
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et al. (1982) found completely analytical solutions for these self-similar configurations. The 
theory applies in the non-relativistic limit when mechanics and gravitation can be approached 
by Newtonian concepts. Written in spherical polar coordinates (r, ^,0), the density profile 
of a Hayashi-Toomre model has the form: 
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where V is the rotation velocity and a is the isothermal acoustic (or stellar dispersive) speed. 
Notice that, apart from a trivial scahng relative to a^, the Hayashi-Toomre models form a 
linear sequence characterized by the single value V^/a^. Notice also that p = for ^ = 
or TT, whereas p — > oo as r — > and the enclosed mass within a sphere of radius r goes 
to infinity as r — > oo. These properties account for our assignment of the name "singular 
isothermal toroid" (SIT) for this generic class of models (see also Li & Shu 1996). 

What is the relativistic generalization of the Hayashi-Toomre sequence when a and V are 
not small compared to the speed of light c? In this paper we confirm the intuitive expectation 
that the linear sequence will broaden into a two-dimensional surface characterized by the 
pair of dimensionless parameters v = V/c and 7 = a^/c^. In the limit = V'^/a'^ ^ 1, we 
further anticipate that relativistic SITs will become highly flattened, like their non-relativistic 
analogs in the same limit; i.e., SITs will become SIDs (singular isothermal disks). ^ 

Relativistic SIDs have been studied by Cai & Shu (2002, hereafter CS), who took their 
inspiration from the cold-disk work of Lynden-Bell & Pineault (1978a,b). CS adopted the 
simplification of negligible disk thickness by assuming an anisotropic pressure tensor, with 
zero and nonzero effective sound speeds in the vertical and horizontal directions respectively. 
Part of the rationale for the present investigation is to ascertain the range of validity of the 
disk approximation when the same effective sound speed a applies in all three principal- 
axis directions, and when the combination v'^/^ is not necessarily large compared to unity. 
In particular, we wish to examine how frame-dragging may distort expectations of toroid 
flatness formed by naive extrapolations from the Newtonian experience. 

The astronomical motivation for the work of Hayashi et al. (1982) and Toomre (1982) 
comes from interstellar gas clouds and disk galaxies. It is unknown whether such objects 
have their relativistic counterparts in astrophysics. For example, active galactic nuclei and 
quasi-stellar objects are increasingly believed to be powered by supermassive blackholes 



"'^In this paper, we use the terminology "isothermal" somewhat loosely to mean that the pressure is 
directly proportional to mass density p, with a constant of proportionality that we call = 70^. Such a 
proportionality can arise, of course, in physical systems in a wider context than when the thermodynamic 
temperature is a strict constant. 
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(SMBHs). Are such SMBHs the central products of slow growth through an accretion disk? 
Or do they result from successive mergers of smaller units (stars or smaller BHs) in the 
nuclei of galaxies? Or could SMBHs have formed from the monolithic gravitational collapse 
of relativistically compact gaseous or stellar systems that resemble SITs? 

Haehnelt & Kauffmann (2001) discuss the difficulties associated with the first two scenar- 
ios. If the last scenario is a viable or even likely possibility, it would be important to develop 
theoretical predictions and observational tests. Such tests provide another motivation for 
the present line of research. In particular, it is known that slowly rotating, non-relativistic, 
gaseous SITs are unstable to inside-out gravitational collapse to form a steadily growing 
pointlike mass at the center (see, e.g., Shu 1977 for the simplest case of the collapse of a 
singular isothermal sphere or SIS). The analogous collapse of a suitably flattened, relativistic 
SIT (during the earhest epochs of galaxy formation) to form a steadily growing SMBH (with 
a Kerr-like geometry) should be accompanied by the copious generation of gravitational 
radiation. Such radiation might be detectable by future gravitational-wave observatories. 

Differentially rotating objects exemplified by SITs and SIDs have other advantages over 
the uniformly rotating, relativistic disks investigated in the pioneering work of Bardeen & 
Wagoner (1971). Prom numerous studies of self-gravitating systems in the Newtonian regime 
(sec, e.g., Toomre 1977, Binney & Tremaine 1978, Lowe et al. 1994, Bcrtin et al. 1989ab, 
Goodman & Evans 1999, Shu et al. 2000), it is known that the transfer of angular momen- 
tum from inside to outside through spiral and barlike density waves is a natural phenomenon 
in configurations that arc differentially rather than uniformly rotating. Moreover, for sim- 
ilar distributions of specific angular momentum, differentially rotating systems have lower 
specific energy, making them more likely states to be encountered in many natural settings 
(Mestel 1963). These considerations add yet more motivation for the continued theoretical 
investigation of relativistic SITs and SIDs. 

The rest of this paper is organized as follows. In §2, we develop the mathematical 
equations governing SITs. In §3.1, we solve the Einstein equations in the Newtonian limit and 
show that the solution thus obtained is identical to those of Toomre (1982) and Hayashi et al. 
(1982). In §3.2 we derive analytical solutions for relativistic SISs when spherical symmetry 
holds. §4 describes the numerical strategy we use to solve the fully nonlinear Einstein 
equations in the presence of finite rotation. §5 summarizes the results of our calculations, 
including the three dimensional velocity field and the degree of fiattening due to rotation. 
Finally, we offer conclusions and commentary in §6. 
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2. Basic Equations 
2.1. Metric 

Unlike the razor-thin disks studied in CS, we wish to consider pressure tensors in this 
paper that are isotropic in three dimensions. In spacetime, the stress-energy tensor is then a 
smooth function of polar angle 6, being positive-definite from pole to equator, with reflection 

symmetry about the latter. Because the sources of gravity are smoother than a delta function, 
the functions, N{9), Q{^)j ^{^) metric adopted by CS, 

ds'' = -rVdt^ + r2e2^-^(# - r^-^e''-''Qdtf + e^-^{dr'' + r'^dO''), (2-1) 

have no kinks as we cross the midplane, unlike the disk case. For clarity, we recapitulate a 
few important properties of this metric. 

In the above equation n is a number between and 1, which we call the gravitational 
index. In our choice of coordinates (but not in Lynden-Bell & Pineault 1978a,b), the radial 
coordinate r has the usual dimensions of length. The presence of the extra factor r" in front 
of each power of dt would then seem to indicate that t does not have the unit of length 
(in geometric units where c = G = 1). This impression is mistaken. As explained in CS, 
the gravitational potential of self-similar SISs, SIDs, or SITs have logarithmic dependences 
on the radial coordinate r. For example, except for slight notational modifications, Toomre 
(1982) and Hayashi et al. (1982) found that the gravitational potential associated with the 
density distribution (1-1) is given in conventional units by 

$ = (V^2 ^ 2a^) Inr + lc^N{e), 

where N{9), up to an additive constant, has a functional form that can be obtained by 
integration of the expression for N' to be given in §3.1. When exponentiated, as in — 
_g2$/c2 _ _^2n^N T^ii]^ ^ _ ^2 _|_ ^j^jg logarithmic dependence yields one troublesome 

factor of r" for each index of t involved in a metric coefficient associated with the differential 
of time. However, an arbitrary length scale L introduced to make the argument of the loga- 
rithmic term formally dimensionless, i.e., to convert Inr to ln(r/L) would add an arbitrary 
constant to the potential $ that would have no physical conseqences. The factors of r" or 
r^" then become the dimensionless combinations (r/L)" or (r/L)^"^, which would restore to 
time t its usual unit. The arbitrariness of L then merely reflects the freedom to scale in a 
problem that is self-similar and lacks intrinsic length scales. For notational compactness, 
we have adopted the convention of setting L equal to unity. In any contour-level figure, 
therefore, the reader can make whatever length assignment she or he might like to any given 
contour, and scale all other contours accordingly. (Except for size, they all look the same.) 
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The coordinate 6 differs from the normal definition of colatitude by a constant factor 
k > 1 (sec again the discussion of CS). If 9 were defined as usual to have a range from to 
TT (as in Lyndcn-Bell & Pineault 1978a,b), k would appear elsewhere in the problem as an 
eigenvalue to be found by satisfying certain boundary conditions. Such eigenvalue searches 
in nonlinear differential equations are numerically very expensive. By absorbing the factor k 
into the definition of wc change the nature of the eigenvalue search to placing the midplane 
in a proper location 9q ^ -k/2. Determining the value of 6q turns out to be possible via a 
computationally non-taxing shooting scheme in our formulation of the numerical solution of 
the Einstein equations. 

Physically, the point is as follows. The gravitation of a flattened distribution of matter 
curves space in such a way as to warp the polar angle 9 if we choose el'^^^^^^^^^'/^r to be 
the radial distance from the origin along any path of constant (f) and 9. (The assumption 
of axial symmetry allows us to preserve the usual meaning of the azimuthal angle 0.) This 
warping is such as to make the total angle coverage from north rotational pole to south 
rotational pole greater than tt (but always less than 27r). 

Apart from the transformations described above, the metric (2-1) has the only form 
consistent with the requirements of self-similarity, i.e., it satisfies the scaling law ds'^ — > a^ds^ 
when r — > ar, and t a^~"-t for a constant a. As explained in CS, the scaling of t arises 
from the gravitational redshift associated with climbing out of the potential of this problem. 
Limiting the effect to redshifts then imposes the allowable range < n < 1, as we noted 
earlier and as will be confirmed in detail later. 



2.2. Einstein's Equations 

We denote ordinary differentiation by 9 with a subscript in this variable, and define 

7l(^a){b) = ''"^e^~^-R(a)(6)! whcrc R{a){h) IS the usual Ricci tensor. Wc adopt as a convention 
that numerical indices or indices in parentheses are tetrad indices, which are raised and 
lowered with the Minkowski metric, while Greek indices are vector indices and arc raised 
and lowered with the metric coefficients of equation (2-1). In this notation, the scaled Ricci 
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tensor resulting from the metric (2-1) has the nontrivial components: 

27^(o)(o) = Ne^e + NgPe + 2n(l + n) - Q^{[{\nQ)e - Pe + Nef + {I - nf] , 
27^(o)(l) = Qe,e + QePe - Q [Pe,e - Ne,e + {Pe - Ngf + Pe{Pe - Ng) + 2(1 - n)] , 
T^mo) - 7^(l)(l) = Pe,0 + Pe" + {n + l)\ 

27^(2)(2) = Ne,e - Ze,e + Pe{Ne - Ze) + 2n(l - n) + Q\l - n)\ 
27^(2)(3) = {ri+ l)Ze - 2niV, + Q\l - n)[Pe - (InQ), - iV,], 

2 [7e(3)(3) + 7e(o)(o) - 7^(l)(l) - 7e(2)(2)] = 2Pe^e - A^e' + 4/1^ + g^jfTv^ + Q^g)^ _ p,]2 _ _ 1)2 

(2-2) 

in the locally nonrotating observer (LNRO) frame defined as usual by 

e(,; = (0,r-^e^^-^,0,0), 

e./ = (0,0,e^(^-^),0), ^^"^^ 



'(2) 
13) 



(0,0,0,r-ie^(^-^)). 



For the matter part, we adopt an isotropic pressure. The pressure p and energy density 
e are related by 

P = le, (2-4) 

where is the sound speed. (In conventional units, e = pc^, and = 7C^-) In 
equilibrium solution, there is no vertical or radial velocity, and we may write in the LNRO 
frame 



Unlike the Newtonian case, where the Poincare-Wavrc theorem ensures that the velocity field 
is independent of z in cylindrical coordinates (e.g., Tassoul 1978), we will have to assume 
here the velocity is a function of the polar angle 9. For a perfect fluid, the stress energy 
tensor is given by 

T{a)(b) = (£ + P)U{a)U{b) + Fn{a){b) (2-6) 

Explicitly, the non-zero components are 



hm) = — ^(o)(i) - — -5-' 

2 (2-7) 

7 + 

" ^ l_^2 ' ^(2){2) = 7£ = 7^(3)(3)- 
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The equation of motion T 



(a)(6) 



m 



— has two nontrivial components: 



(f) : -n + v'^ + 27- 



1 + 7 



+ Qv{l-n)^0, 



{§) : v{Q + v)Pe - vQe - {\ns)e^^—^ - ^7V,(1 + v'' + 2Qv) 

1 + 7 2. 



(2-8) 



0. 



The first equation yields radial force balance: with —n representing gravity; f ^ representing 
centripetal acceleration (a factor 1/r has been cancelled from all terms); 27(1 — f^)/(l + 
7) representing the pressure gradient (the factor 2 comes from differentiating a pressure 
that is proportional to while the factor (1 — v'^)/{l + 7) comes from making various 
inertial corrections); and Qv{l — n) representing the effect of the dragging of inertial frames. 
The second equation describes a similar balance in the 9 direction. In the problem of an 
infinitesimally thin disk, both equations are evaluated on the equatorial plane only. If we 
then impose symmetry about the equator, the second equation is identically satisfied upon 
integrating across the mid-plane. For the 3-D configuration we are considering here, the 
second equation provides a non-trivial consistency relationship for the energy density. 

To proceed further, let us define 



.,Z-N 



e = (1 + n)^, 



(2-9) 



and denote differentiation with respect to © by primes. The Einstein equations may now be 
written as 



Written out in full, the components are 



T 



(2-10) 
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= N" + N'P' + 



2n 
1 + n 



Q" + Q'P' - Q 



g2|[(lng)'-P' + iVf +(i-^^ 

P" _ N" + (P' - N'f + P'(P' - N') + ~ ""l 

{l + ny_ 



2-fe ^ P" + P" + 1, 

(1 - 7)£ = N" - Z" + P'{N' - Z') + ^j;^^ + 



2n 



n 



1 + n 



1 + n 



(l + n)2 
[P' - {\nQ)' - N'l 



n 



1 + n 



47£ = 2P'Z' - N'^ + 



-8- 



With some algebra, we may separate out a part of Einstein's equations and the equation 
of motion (2-8) to put them into a set of "dynamical" equations involving second- and first- 
order derivatives, respectively, of metric and matter variables: 



Q" = -e{2v + Q + Qv^) 



P" = 2j£-l- P'^ 



2n 
1 + n 



n 



1 + 7 



- Q'P' + Q 



1 + n 

1 — n 
1 + n 



(2-lla) 



+ (P' - N'Y - F 



(Ine)' 



v'^P' -vF- -N'{1 + v^] 



1 + 7 



1 — n 



7(1 — v^) 1 + n 



{TV' - QF}, 



(2-llb) 
(2-llc) 

(2-lld) 



where F = QN' + Q' — QP'. This set of ordinary differential equations (ODEs) is supple- 
mented by a set of constraint equations from the rest of Einstein equations and equations of 
motion involving only first- and zeroth-order derivatives, respectively, of metric and matter 
variables: 



0^ + + 
47£ 



1 + n 



1 + 7 



+ Qv{l-n), 



+ P'N' 



1 + n 



+ F^ + 2FP'Q 



n 



1 + n 



n 



1 + n 



(2-12a) 
(2-12b) 



Notice that Z decouples from the other three metric coefficients. Its dynamical ODE 
can be replaced by one of the equations of motion associated with the contracted Bianchi 
identity. The function Z can be easily obtained through integration once the other three 
metric coefficients N, P, and Q are found. 

Equation (2-12a) for radial force balance has the following important implication: v 
cannot be a constant, as in the nonrelativistic case, but must vary with if the coefficient 
governing frame dragging Q does. Moreover, because Q{Q) must vanish on the rotational 
poles, © = and 2©0) v will generally not be zero at the rotation axis for general n and 
7 (e.g., y 7^ for general Hayashi-Toomre SITs). Indeed, because (1 — n) and Q are both 
positive, V will achieve its greatest value on the rotational pole, where frame dragging cannot 
help centrifugal forces and pressure gradients to balance self-similar gravity. 

As a final comment, notice that although equation (2-12b) can be used to compute 
energy density e on computational grid points, it does not manifestly require a positive 
definite value for e. Fortunately, we know that the other non-trivial contracted Bianchi 
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identity guarantees that equation (2- 12b) is a redundant relationship, given the remaining 
independent equations of the governing set. Thus, we can use equaton (2-lld) to integrate 
for In e (always obtaining a positive value e no matter what the sign of In e) , and then check 
numerically that equation (2-12b) yields a consistent result (which it always docs). Thus, 
the remaining equations (2-11) and (2-12a) form a complete set to determine the spacetime 
geometry. 



2.3. Boundary Conditions 

The equations (2-11) form a set of first-order non-linear ODEs in the five variables 
{A^', Q, Q', P', Ine}. Thus, five boundary conditions will specify a solution. As explained 
below, wc have three boundary conditions each that wc desire to impose at the pole and the 
equator, and therefore we seemingly have an overdetermined situation. In fact, because the 
location of the equator is not given a priori, we have a well-posed problem. 

Let us start with the rotation axis. As usual, one wants a circle in the direction 
specified by constant to have vanishingly small circumference if ^ 0. Thus, we need 
— >■ 0, which requires P to go to — oo at the pole = 0. In order to have a non-singular 
geometry on the axis, the Ricci tensor must remain finite as ^ 0. The boundary conditions 
on the axis that are compatible with the Newtonian limit are therefore 

A^' = : The gravitational field is smooth across the axis. 
Q = : There is no frame dragging on the axis. (2-13) 
P' — -|-oo : Coordinate singularity. 

We noted earlier that the fiuid velocity v is generally non-vanishing on the rotation axis 
(as is true, e.g., in the Hayashi-Toomre models). We anticipate therefore that the infinite 
"centrifugal" effect at the pole will drive away all matter from it and create a cavity there 
(explaining, e.g., why the Hayashi-Toomre configurations are "toroids"). If we suppose that 
a near- vacuum situation applies also to relativistic SITs and approximate e to be vanishingly 
small near the rotation axis, we can show that equation (2-lld) has the general solution, 

P = ln[/s:sin(0 - 0i)], (2-14) 

in a small neighborhood of = 0i where e vanishes, with K and 0i being integration 
constants. The coordinate singularity imposed by P' being infinite on the north pole then 
identifies that 0i — 0; i.e., P' ~ 1/0 for small 0. Notice finally that the divergence of 
P' at the north pole is consistent with the vanishing there of the energy density e on the 



-10- 



left-hand side of equation (2-1 Id). Thus, from a physical point of view, we can now see that 
the boundary conditions Q = and P' = +oo at the axis are not two distinct requirements, 
but only one (once we fix the north pole to be at © = 0). 

We recapitulate. The condition Q = at the pole makes v generally nonzero on the 
axis. A non-vanishing v on the rotation axis centrifugally expels matter from the region 

G = 0. This makes e zero, leading to a logarithmic divergence of P(©) at = that allows 
the automatic satisfaction of P' = +oo at the north pole, once we have located it properly 
relative to the equator. The latter can be assured by a proper technique to determine Qq 
(see below). 

Next, we demand the solution to be symmetric about the midplane. This implies that 
the first derivatives of metric coefficients have to vanish there. Thus, at the midplane, 

N' ^P' ^Q' ^ 0. (2-15) 

Notice from equation (2-lld) that e' vanishes automatically once the other three conditions 
are met. Equations (2-13) and (2-15) are the six boundary conditions referred to at the 
beginning of this subsection. 

3. Special Solutions 

3.1. Newtonian Limit 

Before we describe the numerical procedure and general solution, let us pause to discuss 
the Newtonian limit, where v and 7 are small compared to unity. This limit should yield 
Toomre and Hayashi's result. In the Newtonian limit, gij approaches the metric of Euclidean 
space in spherical coordinates; g^i vanishes; and g^oo ~^ — exp(2$) ~ — 1 -|- 0{v^). In our 
metric, this procedure translates to 

7,^;^£,A^ = 0(r^) < 1, © = ©0 = 5, P^lnsin^, Q = Z = 0. 

With these simplifications, Einstein's field equations reduce to the usual Poisson's relation, 
and the equations of motion recover those of Newtonian fluid dynamics. Instead of (2-lld), 
we will revert to using (2-12b) to determine the energy density. Equations (2-12a), (2-12b), 
and (2-lla) now have the approximate forms: 

= -n -y^ 27 

7£ = n^ + ncot^Ar'- ^iV'^ (3-1) 
N" ^e-N' cotB -In. 



- 11 - 



These equations may be combined to give a single nonlinear first-order ODE for N': 



TV" = ( - - 1 )7V'cot^ - -^N'^ + --2n, 
J J 47 7 



(3-2) 



to be solved subject to the boundary condition A?^' = for ^ = 0. The appropriate solution 
is then 



n 



+ 27, 

2 



n 



7 
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27 2 



(3-3) 



N' = 2nl cot e-csce tanh 



n , U 
— m cot - 
27 2 



The first equation recovers the usual Newtonian simplification for mechanical equilibrium: 
gravity n is balanced by centripetal acceleration and the specific pressure gradient 27, with 
no relativistic corrections for inertia or frame dragging. The second and third equations give, 
respectively, the density distribution (proportional to e) and the self-consistent gravitational 
field (proportional to A^') required to achieve mechanical equilibrium in the 9 direction. 
Apart from notational differences, these results are indeed what Toomre and Hayashi found. 

For the nonrotating configuration, we expect a spherically symmetric solution. In this 
limit, n — 2'y, and 

£ = 47, N' = 0. 

Expressed in dimensional variables, this is the familiar result, p — e/(? — a^/27rG'r^, appro- 
priate for a SIS. 



When there is even a slight amount of rotation, n > 27, and e 
for e < 1, 

e (X cx 6"'/'^, 

and curves of constant density p oc e{9)/r'^ follow 



on the axis. In fact. 



(3-4) 



Thus, except for a very narrow range of angles in the meridional plane near the pole, where 
isodensity contours plunge toward the origin, isodensity contours otherwise look nearly cir- 
cular, r ^ constant, when ^ 1. As j'l'y increases because of a greater importance of 
rotational compared to pressure support, isodensity contours become more flattened toward 
the equatorial plane. Indeed, the figures in Toomre (1982) and Hayashi et al. (1982) show 
that the solutions rapidly approach a disk-like solution as one increases the level of rotation 
to 1)^/27 » 1. We expect the same quahtative behavior to extend into the fully relativistic 
regime. 
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3.2. Relativistic SISs 

In the absence of rotation, the general relativistic solution is a spherically symmetric 
one with Q = and v = 0. Equations (2-12a) and (2-12b) then yield 

Equation (2-1 Id) is now trivially satisfied, whereas the right-hand side of equation (2-lla) 
is identically zero, implying that N' is a constant. In fact, N' has to be if we are to satisfy 
the boundary conditions (2-13) and (2-15). With A^', Q, and v all zero, the equation for 
Z' shows that it is also 0. Without loss of generality, we may then choose the constants 
N and Z themselves to be 0, so that the metric (2-1) involves no extra scale factors in an 
interpretation of r as the radial distance from the origin along any path of constant {t, 9, (p). 
Thus, the only nontrivial metric coefficient is P, which satisfies the ODE (2-1 Ic): 

P" + = -a\ (3-6) 
where the positive constant a < 1 is defined as 

« = V 1 - 2^7 = ^^73:^ ■ (3-7) 

Equation (3-6) may be put into the form 

(e^)" = -a^e^, (3-8) 

and may be integrated, subject to the pole boundary condition = at © = 0, to yield 

e^ = Asm{aQ), (3-9) 

with A a nonzero constant. 

When we apply the equatorial boundary condition, P' = at = Go, we obtain the 
identification cos ckOq = 0, or 

e„ = ^. (3-10) 

The only remaining task is to determine the value of the constant A. We accomplish this 
task by requiring that, in the presence of spherical symmetry, half the circumference of a 
great circle in the meriodional plane, 2r9o — 2rQo/{l+ri) according to the metric (2-1), must 
be equal to half the circumference of a great circle in the equatorial plane, 7rre^^^°^ — Anr. 
This equality generates the identification, 

^ = -FT V = ; < 1- (3-11) 

a{l + n) ^1 + 67 + 72 



-13- 



Thus, the ratio of the circumference of a great circle 2A7rr to its radius r from the origin is 
less than 27i for a relativistic SIS with 7 > 0. This is a familar result known for relativistic, 
self-gravitating, pressure-supported (and therefore spherically symmetric) equilibria, but it 
is usually derived for Schwarzschild (or Oppenheimer-Volkoff) coordinates, where 27rrs is 
defined as the circumference of a great circle, and rs is not the true radial distance from 
the origin. What appears as distortions of angles in our metric (e.g., 60 > 7r/2 is the 
"angle" between pole and equator) transforms as distortions of the radial dimension in the 
Schwarzschild metric. If we were to adopt the Schwarzshild description for angles, the radial 
distance r to the origin would turn out to be a power law of rs, with an exponent different 
from unity. 

As a check that our definitions of 9 and lead to no true distortions of angular rela- 
tionships in the spherically symmetric case, we note that the surface area of a sphere can be 
calculated from the metric (2-1) as 

If we compare this formula to the circumference of a great circle 27rr/a{l + n), we note that 
the ratio of the surface area of a sphere to the square of the circumference of a great circle on 
the surface of that sphere equals I/tt, the same relation as for Euclidean geometry (QED). 
Angles (polar and azimuthal) maintain their usual relations; it is only the radial direction 
that suffers true distortion in a relativistic SIS. This completes our discussion of SITs with 
analytically soluble forms. 



4. Numerical Solution 

From dimensional considerations, the solutions for rotating SITs form a two-parameter 
family. The equation of state is a result of microscopic physics, which is independent of the 
macroscopic spacctimc geometry. Thus, the square of the isothermal sound speed, 7, is a 
natural choice for one of the parameters. The other convenient constant of the problem, n, 
measures the strength of the gravitational field. Ideally, we might have preferred to specify 
in advance the amount of rotation. However, this is inconvenient as is a function of "polar 
angle" ©. We therefore fix on 7 and n to parameterize our solutions, and let this combination 
determine implicitly the amount of rotation on the equator and elsewhere. 

To implement a practical numerical scheme, we solve equation (2- 12a) as a quadratic 
relation for v to obtain 

_ -g(l + 7)(1 - n) + V(l + 7)2(1 - n)^Q^ + 4(1 - 7)(n + 717 ^ 

2(1-7) ■ 
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Using the above relation, we compute v{Q) from Q{Q) at every step for use in equations 
(2-11). Causality demands v < 1 and 7 < 1, and Q must have the same sign as v, which is 
taken to be positive by convention. Thus, the solution space is constrained to satisfy 

< 7 < 1, <n<l. (4-2) 

In fact, (1 + 7)7(1—7) times 71 — 27/(1 + 7) is the value of at the pole. 

For convenience, we start at the midplane. Because the governing ODEs, equations 
(2-11) and (2-12), contain no coefficients that depend explicitly on itself, these equations 
are formally invariant to an overall translation of the value of at the midplane. Thus, a 
simple transformation of variables from to ^ = — 0o means that we do not need to 
guess a value of 0o at the equator, but can obtain it afterwards when we reach the pole, as 
defined by the value of where Q = is reached. Attainment of Q = at the pole leads 
to the automatic satisfaction there of P' — +00 as we discussed earlier. 

The boundary conditions (2-15) define the corresponding values of N', P', and Q' at 
the new midplane label ^ = 0. In principle, we need to guess midplane values of Q and e to 
begin the integration toward the north pole. In practice, we can eliminate the need to guess 
e at the equator, and use equation (2-12b) to obtain this value, given 7 and n. Thus, we are 
left with a one- dimensional shooting task of adjusting a guessed midplane value oi Q = q to 
satisfy the boundary condition A^' = at the pole \t' = — 0o- Since we want the local energy 
density to be positive in equation (2-12b), a useful value of q must reside within the interval: 

2n , , 

0<q< . (4-3) 

1 — n 

A slight difficulty enters to complicate the above program that deserves mention. When 
regarded as a function of q, the relevant zero of A^'(O) turns out to be a double root; i.e., 
dN'{Q — Q)/dq — a,t the sought-after value of q. We therefore numerically compute 
dN'{Q — 0)/dq and refine the grid of possible q values when this quantity becomes small, 
so as not inadvertently to jump over the root of A"'(0 = 0). 



5. Results 

To make meridional cross-sectional plots of isodensity surfaces, we use "cylindrical" co- 
ordinates {t, zu, (f), z) that are the counterparts of the system used by Lynden-Bell & Pineault 
(1978a,b): 

w^r^/''sm{9/k), z ^ r^'^ cos{e /k), (5-1) 
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(5-2) 



TT 

is the factor > 1 for the models of this paper discussed in §2. 



In these coordinates the metric reads 



+ k\w^ + z^-h^-^^idw^ + dz^), 



(5-3) 



with the 9 dependences of the metric coefficients N{9), P{0), Q{9), Z{9) now to be inter- 
preted by making the substitution 



Because of spatial curvature, faithful representations of both lengths and angles (in 
the meridional plane) would require curved sheets of paper or a 3-D embedding diagram. 
Conventional publishing limitations restrict our ability to do full justice to such possible 
graphical representations, and we have compromised by making {zu, z) plots of the isodensity 
contours using the coordinate definitions (5-1). We emphasize, however, that when n is not 
small (i.e., when gravity is not weak), such plots give graphical relationships of isodensity 
contours for the coordinate labels only, and do not accurately represent spatial separations 
or angles. Figure 1 gives such plots for the case 7 = 1/3 (e.g., an ultrarclativistic ideal gas 
whether "isothermal" or not) when n ranges from a minimum value n = 1/2, corresponding 
to no rotation (i.e., a SIS), to a value close to the maximum value n — 1, corresponding to 
the most rapidly rotating SIT possible. 

When we introduce rotation, i has to vanish on the result of the centrifugal 

emptying process described earlier. Near the pole, we assume £ oc 0'^ for A > 0, and adopt 
a series expansion for P', using the "dynamical" equations (2-11): 



where we have used equation (2- 12a) on the axis. Notice that in our constrained solution 
space (4-2), A > 0. For a slight rotation, A ^ 1, the power-law behavior of the energy density 
near the axis causes it to change at a normalized radius from 1 to nearly discontinuously 
as the pole is approached. This behavior is depicted in Figure 2. Expressed differently, 
the parameter A controls the slope of isodensity contours near the pole. As A — > 0+, the 
isodensity curve in the meridional plane switches from being horizontal to being vertical very 
quickly as the curve approaches the origin (see Fig. 1). Although e has to vanish formally 
on the axis. Figure 1 shows isodensity contours essentially to be spherical when A is small. 




(5-4) 
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In general, for positive values of 7 < 1, A ranges from for the minimum allowed 
gravitational index n (and no rotation) to 00 for n = 1, which is the maximum gravitational 
field strength allowed by causality. When 7 = 1, a limiting process is required. If we take 
the limit 7 — 1 first, then we are forced to take n = 1, and we recover the nonrotating 
solution of a self-similar sphere with A = 0. However, if wc fix n = 1 and let 7 — > 1, then 
the result is a maximally rotating SIT with X = 00. This ambiguity is a manifestation of 
the sphere-toroid transition at zero rotational velocity discussed earlier. 

From a purely mathematical point of view, the double-limit ambiguity arises because of 
the non- uniform convergence of solutions towards the singular point n = 7 = 1. After all, 
the Einstein equations formally allow solutions with n > 1. However, these spacetimes are 
acausal and have divergent energy density on the axis, and thus are irrelevant for physical 
considerations. We shall not digress further on this issue. In the following analysis, it is 
understood that we always take 7 — > 1 first. 



5.1. Velocity Field 

With the three-dimensional solutions, we have a chance to study the velocity field as a 
function of polar angle. In a disk of infinitesimal thickness, for each value of 7, the rotation 
speed in the plane of the disk is limited to a maximum value that arises when the frame 
dragging becomes infinite: Q — > 00 accompanied by n — > 1. Empirically, CS (see their Table 
1) found the maximum rotational velocity for a disk to be a function of 7 given approximately 

by 

Vc = . (5-5) 

" 2.294 + I.O9I7 ^ ^ 

The fact that the maximum velocity any disk can have is only about 43.8% of the speed of 

light came as a surprise when Lynden-Bell & Pineault (1978b) first discovered it for a cold 

disk, and it has remained a mystery until now. 

Three-dimensional SITs are more reahstic systems than completely fiattened SIDs (with 
their infinitesimal thickness resulting from an artificial imposition of a highly anisotropic 
pressure tensor). In a SIT, as seen in Fig 3, the velocity field increases towards the axis, 
reaching a significant fraction of the speed of light in most cases. This qualitative result was 
anticipated by our previous discussion. The maximum of the velocity field occurs on the 
axis where Q — and is given by equation (4-1) as 

2 n + n7 — 27 



V 



1-7 

For a fixed value of 7 e [0, 1), I'axis increases from for minimum rotation to 1 (the speed 
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of light) for the ultrarelativistic hmit. The sequence is then terminated while the equatorial 
velocities are still a minor fraction of the speed of light, comparable, in fact to the subluminal 
values given by equation (5-5) for the pure-disk solutions. (For more precise comparisons, 
sec below.) The naive expectation that relativity sets limitations on physical systems only 
because matter acquires velocities approaching the speed of light is indeed met for SITs 
where models terminate when (vanishingly small amounts of) gas near the pole become 
ultra-relativistic. It is only SIDs that present an apparent puzzle because the assumption of 
a vacuum above and below the exact midplane of the disk removed all material tracers of 
this ultra-relativistic motion a priori! 

The restriction of equatorial rotation velocities to modest fractions of the speed of 
light implies that differentially rotating, gaseous configurations with isotropic pressure con- 
tributing a significant fraction of the total support against self-gravity, cannot be too highly 
flattened. Despite this caveat, however, it is remarkable how useful the approximation of a 
disk-like geometry can be to predict certain physical characteristics of more realistic systems. 
In Table 1, we tabulate as a function of 7 the maximum equatorial velocity Vc attainable by 
a SID or SIT when the gravitational field acquires the maximum index n = 1 (reached by 
extrapolation from values of n somewhat smaller than unity). From this tabulation, we see 
that the two values of Vc are in very good agreement when 7 -C 1. This is a natural result 
because the critically rotating SIT is extremely flattened at low levels of pressure support. 
However, the good agreement persists to relativistic values of 7 ~ 0.5, a surprise since Vc/7 
is far from being large compared to unity, as is required formally for the disk approximation 
to hold. Indeed, complete breakdown of the disk approximation does not occur in the Table 
until is significantly sma//er than unity. Recalling that the ratio of pressure to energy 
density (including rest mass) is given by 7 = 1/3 for ultrarelativistic ideal gases, we see that 
Table 1 implies excellent agreement in the gross properties of SITs and SIDs with realistic 
equations of state which are maximally rotating. (The apparent discrepancy at 7 = arises 
from an imperfect extrapolation of SITs to this limit.) This finding relieves some of the wor- 
ries that one might otherwise have about the applicabihty of results derived from a flat disk 
analysis when parameters are used that would give such disks appreciable vertical thickness 
if the pressure were isotropic rather than anisotropic (see the caveats expressed in CS) . 



6. Conclusion and Commenteiry 

With the aid of self-similarity, we have constructed a two parameter family of semi- 
analytic solutions to the Einstein fleld equations, parameterized by the gravitational index 
and an isotropic pressure. The isodensity surfaces are toroids, qualitatively similar to their 
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Newtonian counterparts (Toomre 1982, Hayashi et al. 1982). 

The main difference between relativistic SITs and their non-relativistic counterparts 
is that the former do not satisfy, even approximately, the Poincare-Wavre theorem (see, 
e.g., Tassoul 1978). In the Newtonian hmit, Goldreich & Schubert (1967) have shown that 
equihbria which violate the constraint of isorotation on cylinders are unstable on timescales 
that range from dynamical to secular. It is not known whether such instabilities persist 
into the relativistic regime if the cause of the departures from isorotation on cylinders is the 
dragging of inertial frames. This issue deserves further investigation. 

We have also gained insight into how frame dragging figures prominently into another 
puzzle. Lynden-Bell & Pincault (1978b) found, and CS confirmed, that the rotational ve- 
locities of completely fiat, relativistic disks are limited to be less than 0.438 times the speed 
of light. This limitation appears quite peculiar, but when we examine the behavior of three- 
dimensional toroids, the velocity field turns out, as summarized above, to be a function of 
the polar angle. While the equator is still rotating modestly (subluminally) , the velocity on 
the axis can reach the speed of light when an appropriate combination of parameters holds. 
As a result, matter is pushed away from the axis much more strongly than one might have 
suspected by examining only the ratio of rotational support on the midplane to the pres- 
sure contribution 7. As a byproduct, a disklike approximation still holds for many SITs to 
a surprisingly good degree when one might otherwise have expected a complete breakdown 
of such a simplified description. 

The direct applicabihty of relativistic SITs and SIDs to reahstic astrophysical systems 
is questionable. They have inner and outer singularities (infinite central density, divergent 
enclosed mass at infinity) that lead to spacetime geometries (quasi-naked singularity at the 
origin, nonfiat spacetime asymptotically) which are conventionally rejected as inapt descrip- 
tions of reality. However, we note that the issue of the possibility of naked singularities in the 
universe has not been settled in any definitive form (Penrose 1998), and non-asymptotically 
flat spacetimes are a common feature of many standard cosmologies, including the currently 
favored versions of accelerating universes with positive cosmological constants. Moreover, 
for a problem to be interesting in general relativity, spacetime has to possess curvature. 
But self-similar spacetimes that have power law dependences on some radial coordinate are 
guaranteed to run into trouble at small and large radii, since the interesting curvature in the 
body of the problem has to be carried into the origin and to infinity (where they become 
"objectionable"). 

Our point here is not to argue for the direct application of the self-similar models 
of the current line of research to any known astrophysical object, but to point out their 
possible utility in the study of a restricted but extremely interesting subset of problems 
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in general relativity: the role of frame dragging, nontrivial equilibria with and without 
convenient symmetries, gravitational collapse, the formation of singularities, the efficiency 
of the generation of gravitational radiation, etc. Computational results in this field are hard 
to obtain as it is; we should not shrink from applying a technique - self-similarity - that 
has proven itself very fruitful in application to many subjects other than general relativity, 
simply because its assumption produces a few blemishes in otherwise acceptable physical 
models. 
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by a Graduate Fellowship awarded to MJC and by grant AST-9618491 awarded to FHS. In 
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Fig. 1. — Meridional cross section of constant energy surfaces for 7 = 1/3. The contours 
corresponds to 1, 2, 4, and 8 times some reference value. The allowed values of gravitational 
index range from n = 27/(1 + 7) = 0.5 to n = 1. 
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Fig. 2. — Rescaled energy density £ as a function of polar angle for 7 = 1/3. The minimum 
gravitational index n = 0.5 corresponds to a nonrotating sphere with constant e. As n 
increases, energy is redistributed toward the midplane. 
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Fig. 3. — Linear rotational velocity as a function of polar angle for 7 = 1/3. 

Table 1. Critical Velocities for SID and SIT as a Function of Sound Speed Squared 



7 


SID Vc 


SIT Vc 


0.0 


0.438 


0.439 


0.1 


0.415 


0.414 


0.2 


0.398 


0.395 


0.3 


0.381 


0.375 


0.4 


0.366 


0.347 


0.5 


0.351 


0.313 


0.6 


0.339 


0.276 


0.7 


0.327 


0.235 


0.8 


0.316 


0.190 


0.9 


0.306 


0.132 


1.0 


0.298 


0.050 



